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Electron transfer (ET) between primary electron donors and acceptors is modeled in the photosys- 
tem II reaction center (RC). Our model includes (i) two discrete energy levels associated with donor 
and acceptor, interacting through a dipole-type matrix element and (ii) two continuum manifolds of 
electron energy levels ("sinks"), which interact directly with the donor and acceptor. Namely, two 
discrete energy levels of the donor and acceptor are embedded in their independent sinks through 
the corresponding interaction matrix elements. We also introduce classical (external) noise which 
acts simultaneously on the donor and acceptor (collective interaction). We derive a closed system 
of integro-differential equations which describes the non-Markovian quantum dynamics of the ET. 
A region of parameters is found in which the ET dynamics can be simplified, and described by 
coupled ordinary differential equations. Using these simplified equations, both sharp and flat redox 
potentials are analyzed. We analytically and numerically obtain the characteristic parameters that 
optimize the ET rates and efficiency in this system. 
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I. INTRODUCTION 

In photosynthetic complexes of plants, eukaryotic al- 
gae and cyanobacteria, quanta of light excite chlorophyll 
dipoles in the antennas of the light harvesting complexes. 
These local energy excitations are then transferred to 
the reaction centers (RCs) of photosystem I (PSI) and 
photosystem II (PSII), where the charge separation oc- 
curs. During charge separation an electron jumps from 
the donor to the acceptor and both donor and acceptor 
become charged. These electron jumps continue through 
out the whole chain of the redox potential. The char- 
acteristic time-scales of the electron dynamics vary from 
a few picoseconds to milliseconds. The primary charge 
separation occurs on a very short time-scale, of a few pi- 
coseconds [1-5]. Because this time-scale is so short, even 
the room-temperature fluctuations of the protein envi- 
ronment do not destroy the quantum coherent effects, 
which were recently discovered in these complexes |a-[l3j . 

Generally, to describe the motion of an electron be- 
tween the donor and acceptor, different approaches can 
be used 0, [lj-[21|. The leading approach is based on 
the well-known Markus theory, which takes into account 
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both the strength of donor-acceptor interaction, distance 
and the dynamics of the protein environment |14J . In this 
way, the election transfer (ET) rate and the efficiency of 
the ET can be calculated [15|]. At the same time, many 
efforts have been devoted to various modifications of the 
Markus theory. One of them is related to taking into ac- 
count the sinks which represent the continuum electron 
energy reservoirs analogous to the HOMO and LUMO or- 
bitals in the chemical systems [12|, LL3|, [22J (see also refer- 
ences therein) . These continuum manifolds are similar to 
those described by the Weisskopf-Wigner model [23J, [24J . 
These sinks serve as additions to the thermal reservoirs 
in ET. The principal difference between the protein (vi- 
brational) and sink reservoirs is that the former behave 
as the bosonic (electromagnetic) environments while the 
later provide additional electron quasi-degenerate states. 
These sinks increase the entropy for an electron escaping 
into the sinks. As a result, the sinks modify the form of 
the Gibbs equilibrium distribution. Indeed, even in the 
case of a flat redox potential one can expect the acceptor 
to be populated with high enough efficiency. 

In this paper, we consider ET in a model which con- 
sists of two discrete energy states, donor and acceptor. 
A dipole-type matrix element of a direct donor-acceptor 
interaction is included. The donor and acceptor inter- 
act directly with two independent sinks which are rep- 
resented by two quasi-degenerate (continuum) manifolds 
of the electron energy levels. In addition, classical (ex- 
ternal) noise interacts with the donor and the acceptor. 
By taking into account these effects, we describe the ET 



dynamics in different parameter region. We obtain the 
conditions when the system of integro-differential equa- 
tions, which describes generally the non-Markovian elec- 
tron dynamics, can be reduced to much simpler system 
of ordinary differential equations. We obtain analytically 
and numerically the ET for both sharp and flat redox 
potentials, and for different amplitudes of noise. The ad- 
vantage of our approach is that it is a rigorous one, so all 
approximations are controlled and justified. The simul- 
taneous contributions of the sinks and noise are analyzed, 
allowing us to derive optimal ET rates for various regions 
of parameters. 

Our paper is organized as follows. In Section II, we 
describe our simplified model with a single sink interact- 
ing only with the acceptor level, and present an analytic 
solution for the ET, in the absence of noise. In Section 
III, we analyze analytically and numerically the simul- 
taneous influence of the single acceptor sink and noise. 
We derive a closed system of integro-differential equa- 
tions which describe the ET. We also find the region of 
parameters for which a system of ordinary differential 
equations can be applied. In Section IV, we analyze an- 
alytically and numerically the simultaneous influence of 
two sinks and noise. In Section V, we discuss the ob- 
tained results. In the Appendices some useful formulae 
are presented. 



II. MODEL DESCRIPTION 

Our model consists of two protein cofactors in the RC 
(donor and acceptor) each with a discrete energy level. 
In this Section, we consider, for simplicity, only a single 
sink which interacts with the acceptor level (the acceptor 
is also embedded in a sink). The sink can be considered 
either as an additional, third cofactor, or as a part of the 
acceptor. (See Fig. 1.) The first site, denoted by \d), is 
the electron donor, with the energy level, E d . The sec- 
ond site, | a), is the electron acceptor, with energy level, 
E a . We model the sink by a large number of discrete 
and nearly degenerate energy levels, N a 3> 1 (Fig. [T]). 
Then, the transition to the high density of states (in the 
continuum limit) is done for the sink. Note that in some 
situations, the sinks can be considered as approximations 
to HOMO and LUMO orbitals with finite densities of 
states. 

The Hamiltonian of this system can be written as 

H t = E d \d)(d\+E a \a)(a\ + ^(\d)(a\ + \a)(d\) 
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where Ei are energies of the sink levels, and V da = V/2 
(See Fig. [Q) 

Using the standard Feshbach projection method [2 
[28j . one can show that the dynamics of the donor- 
acceptor (intrinsic) states can be described by the follow- 



V 



dv. 



V 



al 



y 



aN, 



1 



FIG. 1: Schematic of our model consisting of donor and ac- 
ceptor discrete energy levels, with the acceptor coupled to a 
sink reservoir with a nearly continuous spectrum. 



ing Schrddinger equation with an effective non-Hermitian 
Hamiltonian, U = U-iW (we set H = 1)|^|: 



•^r =**(*)' 



(2) 



where 



n = s d \d)(d\+s a \a)(a\ + ^(\d)(a\ + \a)(d\) (3) 

is the dressed donor-acceptor Hamiltonian, and W = 
(r a /2)\a)(a\. Here T a is the rate describing the tunnel- 
ing from the acceptor to the sink. (See Appendix A for 
details.) 

Equivalently, the dynamics of this system can be de- 
scribed by the Liouville equation, 



ip=[H )P ]-i{W,p}, 



(4) 



where p is the density matrix projected on the intrinsic 
states, and {W, p} = Wp + pW. 

The solution of the eigenvalue problem for the effec- 
tive non-Hermitian Hamiltonian, %, yields two complex 
eigenvalues: E\,2 = Ao/2 ± 0/2, where \q = e d + e a — iT 
and O = \JV 2 + (e + iT) 2 . Here we denote e = e d — e a 
and r = r o /2. The eigenvalues coalescence in the so- 
called exceptional point (EP) defined by the equation, 
0(V, e, r) = 0. Since O is a complex function of its pa- 
rameters, we obtain two real equations: 5RO = and 
SO = 0. One can show that these equations are equiva- 
lent to s = and V = T. 

Note that, in contrast to the case of an Hermitian 
Hamiltonian, where the degeneracy is referred to as a 
"conical intersection" (known also as a "diabolic point" 
[30l|). the coalescence of eigenvalues results in different 
eigenvectors. At the EP, the eigenvectors merge, form- 
ing a Jordan block. (See the review [3l|, and references 
therein.) 

Choice of parameters. This model involves various pa- 
rameters, whose values are only partially known. Our 
choice of parameters is based on the data taken for the 
ET through the active pathway in the quinone-type of 



the Photosystem II RC 



(Note that the values of large times, t > 150 ps. However, the situation is changes 



parameters in energy units can be obtained by multiply- 
ing our values by H ss 6.58 x 10 _13 meVs. For example, 
e = 60 ps" 1 « 40meV.) 



A. Tunneling to the sink 



drastically for the fiat redox potential. For example, for 
e = 0, we obtain 
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In this section we discuss the ET to the sink. (For de- 
tails see [29|].) We assume that initially the electron occu- 
pies the upper level (donor), pn(0) = 1 and ^22(0) = 0. 
With these initial conditions, the solution of the Liou- 
ville equation (IA22[) for the diagonal component of the 
density matrix is: 
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where Q = \/V 2 + (i + ir) 2 is the complex Rabi fre- 
quency, cos^ = (e + ir)/fi, and sin (9 = V/fi. 

Setting fi = f^+if^ = Vp + *9i where p = U 2 +e 2 — r 2 
and g = 2er, we obtain 



fir 
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■ ± —j=yp + Vp 2 + Q 2 



: ± —j= \j -p + ^/p 2 + q 2 . 
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Using these results, we obtain for /022(f) the simple ana- 
lytical expression: 



022 (*) = 



y2 e -n 



fif) 



2(fi? 



cosh f^f — cos fiii) 



(9) 



We define the ET efficiency of tunneling to the sink as 

7?(t) = 1 - Tr(p(t)) = f Tr{W, p(r)}dr. (10) 

Jo 

This can be recast as the integrated probability of trap- 
ping the electron in the sink [10|, [33| , 



»j(t) = 2r / p 2 2{r)dT 



(11) 



Inserting p 2 2 (t) into (fTTJ) and performing the integration, 
we find that the ET efficiency is given by 



v(t) =1 - 



-17 
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(12) 



In |29| it was shown that, for the sharp redox potential, 
the ET efficiency is rather slow function of time. For 
instance, for e = 60 ps" 1 and lOps < V < 40 ps -1 , the 
ET efficiency approaches a value close to 1 for relatively 



Using these relations, we find that at the exceptional 
point the ET efficiency behaves as 

V(t) = 1 - e- rt . (14) 

For V <C T ( V 3> r) the asymptotic behavior of 77(f) is 



r)(t) 



i-e- rt , u>r, 
i-e- y2t / 2r , u«r. 



(15) 



Comparing the obtained results with (fTS"]) . we conclude 
that the highest ET rate is obtained for the flat redox 
potential at e = 0, and V > T. This result can be in- 
terpreted by using a model of a single spin dynamics in 
an effective magnetic field. The above parameters are 
defined so that the effective magnetic field is oriented in 
the positive x-direction, and V corresponds to the Rabi 
frequency-the frequency of rotation of the spin around 
the a;-axis. Then, the above chosen conditions provide 
a rapid transition from the donor to the acceptor, with 
subsequent tunneling from the acceptor to the sink. The 
results of numerical simulations of the ET efficiency are 
presented in Fig. [5] One can see that, when e = and 
V > T, the ET efficiency can approach a value close to 1 
for short enough times, ~ 2 ps. 



III. NOISE-ASSISTED ELECTRON TRANSFER 
TO THE SINK 

In the presence of classical noise, the quantum dynam- 
ics of the ET can be described by the following effective 
non-Hermitian Hamiltonian (for details see Appendix A) : 

V. =^e„|n)(7i| + J~]\ mn (t)\m){n\ 

n m,n 

+ ^^|m)(n|-7^|2)(2|, rn,n = l,2, (16) 

where A m „(f) describes the noise. We denote by |1) and 
1 2) the donor and acceptor states in the site represen- 
tation, respectively. The diagonal matrix elements of 
noise, X nn , are responsible for decoherence, and the off- 
diagonal matrix elements, X mn (m ^ 71), lead to the 
relaxation processes. In what follows, we use a spin- 
fluctuator model of noise, modeling the noise by an en- 
semble of fluctuators [3J-[3a| . 

The evolution of the average diagonal components of 
the density matrix is described by the following system 




Here E n (z) denotes the Exponential integral [37], A — 
l/ln(7 c /7 TO ), a 2 = x(0) and j m and ~f c (j m < j c ) in- 
dicate the boundaries of the switching rates in the en- 
semble of random fluctuators. The correlation function 
includes, besides the amplitude, a, two fitting parame- 
ters: 7 m and 7 C . Taking into account available theoret- 
ical and experimental data |38l - l40l |. we have chosen for 
our numerical simulations the parameters 7 m and 7 C as 
follows: 2j m = lO^ps" 1 , 27 c = lps -1 . 

Using the Gaussian approximation, we obtain the fol- 
lowing expression for the kernel (see Appendix B): 



V 2 
K(t -t) = — cos(e(i 



where 
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FIG. 2: (Color online) Time dependence (in ps) of the sites 
population, pn(i) and p22(t), (top) and ET efficiency (bot- 
tom) (r = 5ps~ ). Solid and dashed lines correspond to, 
pn(t) and p22{t), respectively. Blue line: V = lOps" 1 , e = 0; 
red line (EP): V — T — Sps" 1 , e = 0; green line: V = 
lOps" 1 , e = lOps" 1 ; black line: V = 5ps"\ s = 20 ps" 1 . 



of integro-differential equations (see Appendix B): 

!<Pll(*)) =~f K{t,t')(( Pll {t')) (P22(t')))dt>, 



and we denote D — \g± — 52 1- The result of integration 
with the correlation function (|2"0|) is given by [36] 



dt 



(p 22 (t)}= K(t,t>)(( Pll (t'))-(p 22 (t')))dt' 
Jo 

-2T{p 22 (t)}, 



(17) 



(18) 



where the average ( ) is taken over the random process 
describing noise, and the kernel, K(t,t'), is given by 

K{t,t') = e-^- 4 ') ((V 2 i(t)V 12 (t')) + (V 2l {t')V 12 {t))). 

(19) 

In the rest of this paper, we restrict ourselves to consid- 
ering only the diagonal noise effects, assuming that the 
noisy environment is the same for both the donor and 
the acceptor sites (collective noise). Then, one can write 
Ai(t) = <7i£(£) and \ 2 {t) = g2^{t), where 51,2 are the 
interaction constants, and £(£) is a random variable de- 
scribing a stationary noise with the correlation function, 
x(t-t') = mW)), given by {H 

X (t) = a 2 A(E l (2 lm T) - ^(27^)), r = \t- t'\.(20) 



\ 4 7l 2 n 4 7c 2 

r t 1 1 

47^~4t| 4t|~47| 
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In Figs. [3] and |3] we present the results of numerical 
simulations for the tunneling rate, V = lps -1 , and dif- 
ferent parameters, V, e, and the amplitude of noise, Da. 
As one can see from Fig. [3l a low level of the noise 
does not improve the ET efficiency rates. However, if 
the amplitude is sufficiently large, the noise significantly 
accelerates the ET to the sink (Fig. [4]). 

In [29| we show that for a sharp redox potential, noise 
can greatly improve the rate of ET to the sink. Our 
numerical results presented in Fig. 01 demonstrate that 
this is true for any redox potential, if the noise is strong 
enough. In Fig. [S] we present the results of numerical 
simulations at the EP and in presence of noise for the 
tunneling rate Y = 5ps _1 (e — 0, V — Y = Sps" 1 ). As 
can be observed, at the exceptional point (EP) the noise 
decreases the rate of the ET. The behavior of the sys- 
tem in the vicinity of the EP is rather complicated, and 
the ET efficiency is sensitive to the choice of parameters 
(Fig. [6]). For instance, for the flat redox potential with 
e = and an amplitude of noise, Da — 60 ps" 1 , the ET 
efficiency approaches a value close to 1 for short enough 
time, ~ 2 ps. 





FIG. 3: (Color online) Time dependence (ps) of the ET ef- 
ficiency (V — F — lps -1 , Da = 5ps _1 ). Blue line: e — 
(EP), red line: e = 20ps _1 , green line: e = 40ps _1 , black 
line: e = 60 ps - . 



FIG. 5: (Color online) Time dependence (in ps) of the ET 
efficiency at the EP in the presence of noise (e = 0, V = T = 
5ps _1 ). Blue line: Da = 0, red line: Da — 20ps _1 , green 
line: Da = 40 ps -1 , black line: Da = 60 ps -1 . 
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FIG. 4: (Color online) Time dependence (in ps) of the ET 
efficiency (V — 10 ps -1 , T — lps -1 ). Solid lines correspond 
to r)(t) in the presence of noise with the amplitude Da — 
40ps _1 . Blue line: e = 0, red line: e = 20ps _1 , green line: 
e = 40 ps -1 , black line: e = 60 ps -1 . Dashed lines correspond 
to r](t), in the absence of noise. 




FIG. 6: (Color online) Time dependence (in ps) of the ET 
efficiency in the vicinity of the EP in the presence of noise 
(e = 0,D(J = 60ps"\r = 5PS" 1 ). Blue line: V = 2.5 ps" 1 , 
red line: V = 5ps _1 (EP), green line: V = 10ps _1 , black 
line: V = 20ps _1 . 



can be approximated by the following system of the or- 
dinary differential equations: 



A. Electron transport approximated by differential 
equations 



Under some conditions, the system of integro- 
diffcrcntial equations, 



j t (Pii(t))=-J Q K(t,t')(( Pll (t'))-(p w (t')))dt', 



dt 



(p 22 (t)) = K(t,t>)((p n (t'))-(p 22 (t')))dt' 
Jo 

-2T(p 22 (t)), 



(24) 



(25) 



dt 

d 

dt 



(p u (*)>=-9t(*)((pil(<)>-(P22(t)», 



(26) 



(p 22 (t)) =9t(«)((pii(<)) - ( P 22(t))) - 2T(p 22 (t)), 



(27) 



where JK(t) = /„ K{r)dT. 

Below we obtain the conditions for which the exact 
system of integro-differential equations (J24I) - (l25l can 
be approximated by Eqs. (|26|) - (|2?p. In the first order 
of the scries expansion, we can write 



( Pll (t>))*(p 11 (t))--( Pll (t))(t-t>), (28) 



d 
di' 



{p 22 {t'))^{p22{t))--{p 22 it)){t-t') 1 (29) 



where, in the same order, the derivatives on the r.h.s. are 
taken from Eqs. pjj) . (1271). 

Using these results, after some transformations we find 
that Eqs. flM]) - flU]) become 

j t {p n (t)) = »(*)«Pii(*)> - <P2a(*)»(l - SHi(t)) 

+ 2rSti(t){p 22 (t)) J (30) 

|(P 22 («)) =»(t)«Pll(*)> - (P22(«)»(l - Wl(*)) 



2r( P22 (i))(i-^ 1 (i)), 



(31) 



where fHi(i) = / tK{t)<1t. From here it follows that the 
system of integro-differential equations (|M| - ([23)1 can be 
approximated by the system of the first order ordinary 
differential equations in the interval of time < t < oo, 
if|/;rif(rH«l. 

Assuming that the correlation function, xW? i s a 
rapidly decreasing function, as t — > oo, we can approxi- 
mate 

" (32i 
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FIG. 7: (Color online) Dependence of 8 on e and Da (T — 
lps" 1 ). 



Performing the integration with the kernel 

K(t) = — cos(er) exp ( - Tr - - °' r 2 

we obtain the following estimate: 

V 2 8 V 2 \8\ t 



(33) 
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where 



S = -l 



— -ex | — 
4p \ 4p 2 



4p 
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exp [ -j—^ |erfc 



erfcf — - 
\2p 



(35) 



where p = Da/\/2, q = T + ie, q = T — ze, and erfc(z) 
denotes the complementary error functiori[37]. 

Using the properties of the erfc function, one can show 
that \8\ < 1 for any choice of parameters T,e and the 
amplitude of noise, Da (Fig. [7J. Consequently, the con- 
dition of validaty of the approximation (l2"rJ)) - (|2"7[) can 
be written as V < Da. This rough estimate can be 
improved greatly for the high level of noise leading to 
V <Da. 

In Figs. [51 IHl we present the results of our numerical 
simulations for T — lps" 1 . As one can see, for V < Da 
there is good agreement between the solutions obtained 
from the system of integro-differential Eqs. (p~T|) - (| 1 8 [) and 
the approximate system of differential Eqs. (|2"6")) - (l2~T)) . 




FIG. 8: (Color online) Time dependence (in ps) of the ET effi- 
ciency, r/(t). The solutions of the system of integro-differential 
equations are presented by solid lines (r = lps -1 ). Dashed 
lines correspond to the solutions of the approximate sys- 
tem of differential equations. Blue line: e = 0, red line: 
e = 20ps -1 , green line: e — 40 ps -1 , black line: e = 60ps _1 . 
Top: V = 30ps"\ Da = Sps" 1 . Bottom: V = lOps" 1 , 
Da = 40ps _1 . 

The advantage of using differential equations instead 
of integro-differential equations is (i) the ability to intro- 
duce effective rates and (ii) the ability to compare the 
results with the Markus theory. The computation of the 
asymptotic ET, $H r = limf->.oo «H(t), yields J29| 



m r =- 



UV27T 

8Da 



exp 



exp 



( r+ ^ 2 W r+i£ 



2D 2 a 2 



\V2Da 



2DW J \ 



^2Da, 



(36) 




and e, the rate 9V reaches its maximum value when the 
amplitude of noise Da f=a e. 

Inserting $H r instead of VK(t) into Eqs. (|26f and |27| . 
we obtain the following system of the differential equa- 
tions: 



d 

Jt 

j t (p 22 (t)) =Wr«Pii(t)> - (P22(t))) - 2T( P22 (t)). 



-( Pll (i)) = - 9t r ((pii(*)) - (fla(t)», (37) 

2r(p 22 (t)) 

(38) 



FIG. 9: (Color online) Time dependence (in ps) of the ET 
efficiency. Solid lines correspond to r){i) obtained from the 
system of integro-differential Eqs. (|17p ~(|18 p and dashed lines 
present the results for r/(t) obtained from the system of differ- 
ential Eqs. pg ]l -p7 |l . Blue line: e = 0, red line: e = 20ps~\ 
green line: £ — 40 ps -1 , black line: e = 60 ps -1 (r = lps -1 , 
V — Da = 40 ps -1 ). One can observe the excellent agreement 
between the solutions. 
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FIG. 10: (Color online) Blue line describes the time depen- 
dence of the ET rate, 9i(t). Red line corresponds to the 
asymptotic rate, 9^r, given by Eq. (|36[). From top to bot- 
tom: e = 0, 20, 40, 60 ps -1 (r = lps -1 , V = 10 ps -1 and 
Da = 40ps -1 ). 

In Fig. 1101 we compare the results of numerical calcu- 
lations of the relaxation rate, $H(£) (blue line), with the 
asymptotic expression JHp (red line) for different choices 
of parameters. As can see in Fig. [Til for the given V 
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FIG. 11: (Color online) Dependence of the asymptotic rate, 
£Kr, on e and Da (V = lOps" 1 , T = lps" 1 ). 



The solution is given by 
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where 9^,2 = ^r + T ± \/ml + T 2 . The computation of 
the ET efficiency yields [29} 



v(t) = i - 



(sfti+m 2 )* 



cosh 



(m!-m 2 )t Wi 



5Ki-$H 2 



<H 2 . , (Wi - 9t 2 )t 

■ sinn 



(41) 



r 



Its asymptotic behavior is 



T)(t) 



9*1 



-SRat 



<Ki-<k 2 



(42) 



As one can see from Eq. (J4T|) . there are two ET rates, 
?ti and $K 2 - However, the asymptotic behavior of the ET 
efficiency is defined by the lowest ET rate, 9\ 2 . 

As shown in Fig. [TU1 the ET rate reaches its asymp- 



totic value, D\(t) —> D\r, quite rapidly, at t « 0.1 ps. 
This allows us to use the analytical solutions to describe 
tunneling to the sink with very high degree of accuracy. 
This conclusion is confirmed by our numerical simula- 
tions presented in Fig. [12j One can observe the excellent 
agreement between the ET efficiency given by formula 
flUJ) and the results obtained from Eqs. fl25J) - (|2"7|) . 




The Hamiltonian of the system can be written as 
H t = E d \d)(d\ + E a \a)(a\ + ^(\d)(a\ + \a)(d\) 

+ J2(Ei\i){i\+V di \d){i\+V id \i){d\) 



i=l 



+ Y,(E j \j){j\ + V aj \a){j\ + V ja \j){a\), (43) 
i=i 

where Ei (Ej ) are energy levels of the sinks coupled with 
the donor (acceptor), and V da = V/2. 

After the transition to the continuum spectra of the 
sinks, the system is governed by the effective non- 
Hermitian Hamiltoniancan, % = % — iW, where 

H= £l |l)(l|+e 2 |2)(2| + |(|l)(2| + |2)(l|), (44) 
where H is the dressed donor-acceptor Hamiltonian and 



w=i(r 1 |i)(i| + r 2 |2)(2|). 



(45) 



Passing from Eq. (g3) to Eqs. flU) - (gSJ) we have 
changed \d) —> |1) and \a) -* |2) . (See Appendix A.) 



FIG. 12: (Color online) Time dependence (in ps) of the ET 
efficiency, r){t) {Y = lps" 1 , V = lOps" 1 , Da = 40ps _1 ). 
Top. Left: analytical solution (dashed lines). Right: results of 
numerical simulations (solid lines). Blue line: e = 0, red line: 
e = 20 ps -1 , green line: e = 40 ps -1 , black line: e = 60 ps -1 . 
Bottom. The both graphics are overlapping. 



V 



da 



V 



dl 



V 



dN, 



V 



aN a 



1 

N d 

1 
N a 



IV. MODIFIED MODEL WITH TWO SINKS 



In this Section, we generalize our model by including 
two sinks interacting independently with the donor and 
acceptor. (See Fig. 13.) The main reasons for this 
generalization are the following. When using a single 
sink which interacts with the acceptor, the asymptotic 
for large times of the ET efficiency, rj(t), approaches the 
unit independently of the parameters of the system, and 
the tunneling rate, T. The additional sink, which inter- 
acts with the donor, describes the leakage of the electron 
from the RC in the process of the ET. By manipulat- 
ing the densities of the donor and acceptor states, or the 
tunneling rates, F\ and T2, (see below), one can fit the 
asymptotic behavior of the ET efficiency, 772 , on the ac- 
ceptor (0 < 772 < 1) with its experimental value. This 
approach is used for modeling the dynamics of the ET 
in photosynthetic complexes. (See, for example, |12l[l3l|. 
and references therein.) 



FIG. 13: Schematic of our modified model consisting of donor 
and acceptor discrete energy levels, with the donor and accep- 
tor coupled to independent sink reservoirs with nearly contin- 
uous spectrum. 

The dynamics of the system is described by the Liou- 
ville equation, 



ip={H, P }-i{W,p}, 



(46) 



Further, we assume that initially the electron occupies 
the upper level (donor), pn(0) = 1 and ^22(0) = 0. With 
these initial conditions, the solution of the Eq. (j46| for 
the diagonal component of the density matrix is given by 



Pn(t) =e 



-rt 



Sit . „ . fit 

cos 1 cos sin — 

2 2 



022 (*) =' 



,-r* 



Qt 



sin a sm ■ 



(47) 
(48) 



where T = (ri + r 2 )/2, ^ = \/V 2 + (e + iT) 2 being We define the ET efficiency of trapping the electron in 

the complex Rabi frequency, cos9 = (e + iT)/il, and the acceptor's sink as 
sin0 = V/fi. 

Setting 51 = ill + zf2 2 , we obtain for P22(t) the simple /"* 

analytical expression: " W) =T 2 J pii(r)dT. (50) 



P22W 



^e 



2„-rt 



. ^ cosh fi 2 i- cos fiii). (49) Inserting /9 22 (*) into ^SDJ) and performing the integration, 

2(S2 1 + S2 2 J we obtain 



m{t) 



ri + r 2 



-17 



r(ni + n 2 



-((r 2 + fii)(rcoshfi 2 t + si 2 sinhrj 2 t) - (r 2 -£7 2 )(rcosftit-f!ismrM)) ). (51) 

1 



From here it follows ry 2 (t) — > 770 , as t — > 00, where 

r 2 



'/o 



ri + r 2 



(52) 



A. Noise-assisted electron transfer 

In the presence of classical diagonal noise, described 
by A„ = g n £(t), the effective non-Hermitian Hamiltonian 
can be written as 

K=J2( £ "+9nm-^)\n)(n\ + ^ J2 HH- 



m^n 



(53) 



We describe the evolution of the system by the follow- 
ing differential equations: 



d_ 
dt 
d 



;(p n ) = - m r {( P n) - (p 22 » - ri<p u ), (54) 

di (p 22 ) =9tr«Pn>-<p 22 ))-r a (p 22 ), (55) 

where £Hr = limt->oo 9t(t) is given by 



The motivation to use this simplified description is as 
follows. As was shown in Sec. Ill, approximation of 
intcgro-diffcrcntial equations by the system of differential 
equations is valid for v < Da. In addition, since the 
ET rate reaches its asymptotic value, IH(i) — > 9tr, quite 
rapidly, at t sa 0.1 ps, we can use the asymptotic rates, 
£Hr, instead of 5H(t). This allows us to use the analytical 
solutions to describe tunneling to the sinks with very high 
degree of accuracy. 

The solution of Eqs. ([54"]) - (f55"j). with the initial con- 
ditions, (pn(0)) = I and (p 22 (0)) = 0, is 



(PnW) 



(p 22 (*)> 



2 2^H 2 + A 2 



A 



-•Hit 



2 2^<K 2 + A 2 



-«2* 



SKr 



2 V«^ + A 2 



e -«2t _ g-«l* 



(57) 
(58) 



"* r 8Da [ CXP { 2D 2 * 2 ) CTtC {V2Da 

\T-1e) 2 



exp 



2D 2 a 2 



crfc 



\V2DaJ 



(56) 



where 9ti, 2 = 5H r + r ± ^/mf+A 2 and A = (r 2 -ri)/2. 



The computation of the ET efficiency of tunneling in 
the acceptor's sink yields 



r 2 «H r L _ (gti+gt 2 )* / , 



smh 

2 JRi - 9t 2 2 



(59) 



From here it follows, as t — > 00, 
f?2 (*) -» ??r 



Comparing the obtained results with the ET efficiency 



rVHr 
<Hi<K 2 ' 



(60) 
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without noise, 770 
obtain 

Vr 



■- r 2 /(Ti + T 2 ), (see Eq. ((52 
2Vd\ r 



' 2r$n r + r 2 - a 2 ' 

This can be represented in the form 

(i\ + r 2 )m r 



Vr = VO 



(Tx + r 2 )*n r + r^ ' 



(61) 



(62) 



From here it follows, that when 9t r > r^/^i + T 2 ) 
the ET efficiency ry r « 770- Generally, the ET efficiency 
with the presence of noise cannot exceed the ET efficiency 
without noise, r\ r < 770. 




FIG. 14: (Color online) Time dependence (in ps) of the ET 
efficiency, r]2(t) (Fi = lps -1 , T2 = 5ps _1 ). Blue line: e = 0, 
red line: e = 20 ps -1 , green line: e = 40 ps -1 , black line: 



e = 60 ps" 1 
V= 10 ps" 1 



. Left: V = 40ps~ 
Da = 40ps _1 . 



Da = 40 ps" 



Right: 



In Fig. Qj] we present the results of numerical simula- 
tions of the ET efficiency when two sinks are taken into 
account. As one can see from Eq. (fBTj) . the asymptotic 
value for the ET efficiency, for chosen values of Ti^, must 
satisfy the condition, r] rc < 770 « 83.3%. As one can see, 
the asymptotic value of 77^ is close enough to 770, for cho- 
sen in Fig. [H] (left), parameters corresponding to black 
line. We also would like to mention that, as the results 
presented in Fig. [14] (left) demonstrate, the saturation of 
77 2 (t) happens fast enough, at t w 2 ps. 



V. CONCLUSION AND DISCUSSIONS 

In this paper, we analyzed analytically and numeri- 
cally the simplest model of the electron transfer (ET) 
between two protein sites, donor and acceptor, in the 
presence of classical (external) noise that is character- 
ized by its amplitude and its correlation time. The noise 
is described by the well-known model of two-level fluc- 
tuators. We also include in our model two sinks which 
are represented by quasi-degenerate manifolds of electron 
energy levels. The sinks are directly coupled to the donor 
and acceptor states. This is done using the well-known 
Weisskopf-Wigner model 23, 24]. 

When both noise and sinks influence the ET, the elec- 
tron dynamics becomes rather complicated, and is gen- 
erally described by a system of integro-differential equa- 
tions. We derive the conditions for which a simplified 
system of ordinary differential equations can be used. 



Our approach is rigorous in the sense that all approx- 
imations are controlled and justified. We obtain analyti- 
cally and numerically the optimal ET rates and efficiency 
for both sharp and flat redox potential. Our results are 
useful for analyzing and engineering optimal properties 
of photosynthetic bio-complexes. 
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Appendix A: Non-Hermitian Effective Hamiltonian 

We consider the time-dependent Hamiltonian of a N- 
level system coupled with independent sinks through 
each level: 

JV 

«(*) = X)e»(t)|n)<n| + J2 f3nm(t)\n)(m\ 

n— 1 m^n 

N N n 

+ Y, £(^J^>w+^i n >^i+ l ^»i*»x n i)' 

n— 1 i n — 1 

(Al) 

where m,n — 1,2,..., A. We assume that the sinks are 
sufficiently dense, so that one can perform an integration 
instead of a summation. Then we have, 



H(t) = £e„(t)|n)(n| + £ p mn (t)\n)( 
+ J2( [a n (E)\n){E\g n (E)dE + h. 
+ ]T J ' E\E)(E\g n (E)dE 



where g n (E) is the density of states, and V m 
With the state vector written as 

M = £ (cn(t)\n) + [ c nE {t)\E)g n (E)dE 



the Schrodinger equation, 

. dm)) 

at 



= nm)), 



(A2) 
a n (E). 

(A3) 

(A4) 
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ic nE {t) =Ec E (t) + a n (E)c n (t). 



takes the form 

ic n (t) =E(t)c n (t) + Y, Pn m (t)c m (t) 

/■OO 

+ / a* n (E)c nE (t)g n (E)dE (A5) 

(A6) 

In order to eliminate the continuum amplitudes from 
the equations for the discrete states, we first apply the 
Laplace transformation: 

/>oo 

Cn(t) = I e- st c n (s)ds, (A7) 

(A8) 

Then, from Eq. (IA6|) we obtain 

(s + iE)c nE {s) = -ia n (E)cn(s). (A9) 

This yields c E (s) = — ia„(£ , )c 2 (s)/(s + iE). Inserting 
this expression for c E {s) into Eq. (|A5I) . we obtain the fol- 
lowing system of integro-differential equations, describing 
the non-Markovian dynamics of the TLS, 

ic n (t) =Ecn(t) + Y f3 nm (t)c m (t) 

st, f\a n (E)\ 2 g n (E)dE n 



c nE (t) = e s c nE (s)ds. 

Jo 



i / c n (s)e ds 
Jo 



s + iE 



(A10) 

To proceed further, we change the variable in the last 
integral, s — > —E', so that 



\a n {E)\ 2 g n (E)dE 
s + iE 



-> -% 



\a n (E)\ 2 g n (E)dE n 



E-E> 



(All) 



The next step is to use the identity 



1 



V 



1 



x — x' + iO l r - ,. 

where V = Principal value. This yields 

\a n {E)\ 2 g n (E)dE 



iirS(x - x'), (A12) 



E _ E , = A{E')--T n {E'), (A13) 



where 



A(E') = V I 



l (E)\ 2 g n (E)dE 
E-E' 



(A14) 



T n (E') = 2ir J \a n (E)\ 2 g n (E)S(E - E')dE. (A15) 

Now using the Weisskopf-Wigner pole approximation, 
we evaluate the integrals as follows 23, 24, 41]: 



\a n (E)\ 2 g n {E)dE 



E-e n 



(A16) 



A(E')^A(e n )=V J 

T n {E') » r„(e„) = 2vr / \a n (E)\ 2 g n (E)S(E - e n ) dE 
= 27r 3 „(e n )|a„(e„)| 2 . (A17) 



The Weisskopf-Wigner pole approximation basically cor- 
responds to the assumption that the coupling constant 
to the continuum is a smoothly varying function of the 
energy, e.g. the continuum is treated as a single discrete 
level. 

Inserting (IA16[) into Eq. (|A5[) . we obtain 



ic n (t) = e n (t)c n (t) + ^ Pnm(t)Cm(t) ^c n (i), 

(A18) 



where F„ = T n {E n ) and e{t) = e„(t) - A{E n ). 

Writing \iptf) = X)n c ™Wl n )> we ^ n< ^ * na ^ ^ ne dynam- 
ics of the iV- level system interacting with the continuum 
is described by the Schrodinger equation, 



■ d\i/jN(t)) £, , ,,. 
* -Q t = ri\ipN[t)), 



(A19) 



where % = % — iW is the effective non-Hermitian Hamil- 
tonian, 

H = J2 £ n\n)(n\ + J2 (3mn{t)\m)(n\ (A20) 

being the dressed Hamiltonian, and 

W = lj2T n \n)(n\. (A21) 



Equivalently, the dynamics of this system can be de- 
scribed by the Liouville equation, 

ip = Up-pitf = [H,p]-i{W,p}, (A22) 

where p is the density matrix projected on the intrinsic 
states, and {W, p} = Wp + pW. 

In particular case of the two-level system considered 
in this paper, the effective non-Hermitian Hamiltonian 
takes the form 



u 



1 ( 2e x - iTi V 

V 2e 2 - tT 2 



(A23) 



Comments. The results of this section can be obtained 
using the standard Feshbach projection method [24IJ28J . 



Appendix B: Equation of motion for the average 
density matrix 

In this Appendix, we derive from the Liouville equa- 
tion, ip = [H,p] — i{W,p}, the equation of motion for 
the average density matrix. We will use the interaction 
representation. Considering the off-diagonal elements as 
perturbations, so that % = "Ho + V(t) — iW, where Ho = 
E n Zn\n)(n\+j: n \ nn (t)\n)(n\, V(t) = E m ^n( V mn + 
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Am»(t))|m){n|, and W = (Ti/2)|1><1| + (T 2 /2)|2){2|, we and 
obtain the following equations of motion: 

Pn =*(/5i2V 2 i - V12P21) -Tipn, (Bl) 

P22 = «(/02lVl2 - V2IP12) ~ T 2 p22, (B2) 

P l 2=iV 12 (p ll -p 2 2)-Tp 12 , (B3) V- = r(e i J"o^o(r)«*r )Vr(e -i/ *Ho(r)«fc- )< (B6) 

/5 21 = iV 2 i(pn -P22) -rp 2 i, (B4) 

where r = (rx+r 2 )/2, 

p = T(e^o^M^ )pT(e -/ 'H oMd . ); (B5) Ug . ngEqs ^ _ (E])) we obta . n 

I 



Pxx{t) =/5n(0) + t I e-r^XfaMVtxtf) - V l2 {t')p 2l {t'))dt' , (B7) 

Jo 

p 22 (i) =p 22 (0) + j / e- ra(t -*' ) (/52i(*')^i2(* / ) - ^2i(t')Pi2(i'))^'. (B8) 

Jo 

M*) =pi 3 (0) + t / e- r(4 - t ' ) T>i 2 (t')(-5ii(i') ~ P22(t'))dt' , (B9) 

Jo 

p 21 (t) =P2i(0) + 1 f e-^^i^X/M*') - /5 22 (0)^. (BIO) 



I 

We assume that initially pi 2 (0) = p 2 i(0) = 0. Now, into account that /5n = pn and p 22 = p 22 , we obtain the 
inserting (IB7|) - (|B10|) into Eqs. (|B1[) - (|B4I) . and taking following system of intcgro-differential equations, 



Pn(t) = - / e- r(t - 4,) (v 21 {t)V 12 (t') + V 21 {t')V 12 (t)) (pn(t') - P 22{t'))dt' - T lPll (t), (Bll) 

p 22 (t) = / e-r^- 4 ') (F 2 i(t)Fi 2 (t') + V 21 (t')V 12 {t)) (pu(t') - p^t'^dt' - r 2 p 22 (t), (B12) 

Pia(*)=- / (e-W-V +e- r ^ t - 4 '))(p 21 (t')/5i 2 (i')-^i2(t')P2i(t'))v r i2(t)^'-r Pl2 (t) 





+ *^1 2 (*)(P11(0)-P22(0)), (B13) 



&i(t) = - / (e- r ^ i - t ') + e-^f*- 4 ')) (V 21 (t')p 12 (t') - V 12 (t')P2i(t'))v 21 (t)dt> - T P2l (t) 

+ iV 2 i(t)(pu(0)-p 22 (0)). (B14) 



For the average components of the density matrix this 
yields 
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j t (Pn(t))=-J e- r(t - t ' ) ((^i(t)Fi2(0+^i(t')^i2(*))(pii(t')-P22(*')))* , -ri<pii(t)>, (B15) 

|(P 22 W) = f e- r ( t - t ')((v 21 (t)^i 2 (i') + V 21 (t')V 12 (t)) (pn(i') - P22{t')))dt' - T 2 (p 22 (t)), (B16) 

!<Pi 2 (*)> - - £ (e- ri(t - f,) + e- r ^ i -*'))((F 21 (i')pi 2 (*') ~ ri 2 (OP2i(* , ))^a(*))d« / - I'M*)) 

+ i(^i2(t))(pii(0)-p 22 (0)), (B17) 

|<P2i(*)) = - / (e- ri(i -*'^ + e- r =(*- t '))((F 21 (i')pi 2 (* / ) - ^(f )Pai(*'))%i(*))d« / - r(p 2 i(*)) 

+ i(V2i(t))(pu(0)-p22(0)), (B18) 

I 

where the average ( ) is taken over the random process and so on. Using these relations, we obtain the following 
describing noise. system of intcgro-differential equations for the average 

In the spin-fluctuator model of noise with the number components of the density matrix, 
of fluctuators, TV 3> 1, one has the following relations for 
the splitting of correlations 36] , 

((V 2l (t)V 12 (t') + V 21 (t')V 12 (t))(p n (t') - p 2 2(t'))) = 

((V 21 (t)V 12 (t')) + <V 2 l(t')^12W))«Pll(i')) - (P22(t'))), 

(B19) 



J t (Pn(t)) = " I e- r(t -'' ) «V 2 i WMt')) + {V 21 {t')V l2 {t))) ((p u (f )> - (M*')>K - iMpnW), (B20) 

^{p 22 (t))=J*e- r ( t - t '\(V 21 (t)V 12 (t')) + ^ (B21) 

d 



-{~ Pl2 {t)) =i{V l2 {t)){ Pll {0)- p 22 {Q))- j (e- r ^-^ +e- r ^-^)(VMt)V 21 (t'))(p 12 (t'))dt' 

+ J (e- ri(t -°+e- r2 ( t - t '))(U 12 (t)l/ 12 (t')>(P2i(i')>^'-r(p 12 (i)>, (B22) 

J t (P2i(t)) =i(V2i(t))(Pu(0)- P22 (0))- f (V ri(t -° +er T ^'^){V 21 {t)V 21 {t')){p 12 {t'))dt' 

+ f (e- ri(f -°+e- r2 ( t - t '))(V 21 (t)f 12 (t'))(P2i(i'))^'-r(p 2 i(i)). (B23) 



o 

r 



Using these relations, we obtain the following system nents of the density matrix, 
of integro-diffcrcntial equations for the diagonal compo- 

I 



!<Pii(*)>=-j[ K(t,t')((p 11 (t'))-(p22{t')))dt'-T 1 {p 11 {t)}, (B24) 

j t (p 22 (t)) ^ J\(t,t')({pii(t')) - { P22 (t')))dt' -T 2 (p 22 (t)), (B25) 
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where the kernel, K(t,t'), is given by 

K(t,t , )=e-^ t - t '\(V 21 (t)V 12 (t')) + (V 21 (t')V 12 (t))). 

(B26) 

For the diagonal noise, so that A TO „ = 0, (m 7^ n), 
the kernel can be recast as 

Kit -t') = \ a»(e(t - t'))e-^ t - t '\e lK ^ t ' ) ), (B27) 

where n(t — t') = D J Q £,{t)cIt and D — \gi — g 2 \. 



In the Gaussian approximation, the generating func- 
tional becomes 



/ e iK(t-t')\ =e -(K 2 (t-t'))/2_ 



(B28) 



t-t' 



n 2 (t - t')) = 2D 2 / dr' / dr" X (r' - t"). (B29) 



With help of Eqs. (JB28| - (|B29|) we obtain 



if (i -t') = ^- cos(e(t - t')) exp ( -T(t- t') - D 2 [ dr' f 
2 V Jo Jo 



dr' I dT" X (r' - t") 




(B30) 



Using the results obtained in Sec. Ill, one can show 
that for V < Da the system of intcgro-diffcrcntial equa- 
tions (|B24[) - (|B25|) can be approximated by the following 
system of ordinary differential equations: 



where p = Dcr/y2. q = T + ie, q = T — ie, and erf (z) is 
the error function [371 ] . 



-( Pll > = - W(t)«pii) - (p 22 )) - ^(pn), (B31) 
^(p 22 > =5H(t)((pii>-(p 22 »-r 2 (p 22 ), (B32) 

where 9t(t) = /„ TK{r)dT. Performing the integration 
we obtain 



^U + P t)-crf(. q 



+ ^ exp (^)( erf (i + ^) _crf 



(!))■ 

(B33) 
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